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I. INTRODUCTION 



A topic of general interest is the numerical solution of partial differential equa- 
tions (PDEs). such as the Navier-Stokes equation (Equation II. 1). Many aspects of 
ilie.se types of problems present interesting challenges to constructing numerical so- 
lutions. such as non-linearities, high-speed flows that cause numerical peculiarities, 
unusual boundary conditions, uncommon geometries, and so on. Developing pro- 
cesses that solve these types of problems is difficult, and often results in complicated 
solutions schemes which are costly to implement. A common goal is to attempt to 
sireamline existing solution processes, or develop new ones, in order to reduce com- 
putational complexity and cost. 

At issue is how to transform a continuum problem into a discrete one (dis- 
cretization), and how to construct a numerical process that will solve the discrete 
problem. The finite volume element (FVE) method (see [Ref. 1]) has proven to be 
a useful tool in developing discretizations, particularly for problems that require the 
enforcement of conservation laws. One of the more successful methods for streamlin- 
ing the solution of PDEs, particularly elliptic equations, is a multilevel technique, to 
wit. multigrid (see [Ref. 2], [Ref. 3], [Ref. 1]). While this method enjoys remarkable 
success in solving certain classes of both linear and non-linear PDEs, its applicability 
to solving other types of equations awaits development. The goal of this work is two- 
fold: to apply the FVE method to discretize a particular equation, and to implement 
multigrid in solving the resulting discretized problem. 

The approach that is taken is to begin with a specific physical problem which 
results in the Navier-Stokes equation, and consider a subsidiary problem, namely the 
heat equation; due to the nature of the problem, cylindrical coordinates are used. 
The resulting problem, to which the FVE discretization and multilevel solution are 
applied, is the axisymmetric heat equation, involving the use of a point source. More 
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specifically, finite differences are used to discretize the temporal portion of the prob- 
lem. resulting in the use of a fully implicit backward time-stepping scheme, and the 
FVE method is used to discretize the spatial part. The application of the FVE 
method to a problem in cylindrical coordinates is new. and results in stencils which 
are analyzed extensively. Once the discretization has been constructed, a number of 
relaxation schemes, including both Jacobi and Gauss-Seidel, are considered; a thor- 
ough analysis of these schemes is done, using both the spectral radii of the iteration 
matrices and local mode analysis. The goal is to choose one of the relaxation schemes 
fur use in an iterative solution of the problem; Gauss-Seidel is the method of choice. 
The specifics of the multilevel technique are then developed and analyzed. Local 
mode analysis is performed on the components of the amplification matrix, resulting 
in the two-level convergence factors for various combinations of the operators. The 
multilevel solution process is implemented by using multigrid V-cvcles; the iterative 
and multilevel results are compared and discussed in detail. The computational sav- 
ings resulting from the multilevel process are then discussed. Insofar as multigrid is 
quite successful in a number of different contexts, the results of our work are some- 
what disappointing in that we are not yet able to achieve the same level of success. 
On a more positive note, this work applies the FVE method to a problem in cylindri- 
cal coordinates, and makes use of Maple software to compute the necessary integrals 
to construct the discretization (see Chapter III and Appendix B). Additionally, the 
techniques that are applied do result in a solution to the problem. 
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II. 



PROBLEM STATEMENT 



A. BACKGROUND 



Our ultimate iroal is to be able to numerically solve the Navier-Stokes equation. 



where u is the velocity vector, t is the time, p is the density, p is the pressure, and // 
is the kinematic viscosity. One approach to such solutions is to consider this problem 
d s a collection of smaller problems. By solving each of the smaller problems, we 
can assemble the collection of solution processes to solve the original problem. The 
purpose of this work is to focus on an initial step in the eventual construction of a 
numerical solution to Equation 11.1. 

One way to subdivide a problem such as this into smaller problems is to 
consider a specific physical problem that gives rise to the Navier-Stokes equation, and 
then isolate the different aspects of that problem. To that end, we consider a semi- 
infinite block of metal, with a heat source applied to the horizontal surface; above the 
horizontal surface is an (inviscid) nonconducting gas. The result is a pool of molten 
metal with a free surface, surrounded by the solid portion of the unmelted block (see 
Figure 1 ). The total heat flux Q is constant, and far away the solid approaches the 
uniform cold temperature, T c . The resulting thermal and flow fields are assumed to 
be axisymmetric and steady, and are governed by conservation of energy in the solid 
and bv conservation of energy, momentum, and mass in the pool. We therefore have 
the following system of equations: 




— Vp + 
P 



(II. 1) 



solid : - 7 -- = kV 2 T 

dt 



(II.2) 



liquid : — b u • VT = kV 2 T 



dt 



(II.3) 



du _ 1 2 _ 

— — b u • Vu = — Vu + 

Ot p 



(11.4) 

(11.5) 



V • u = 0 
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where T is the temperature, k is the thermal diffusivitv, and u. t , p , /;, and u are as 
above. 

Q 




Figure 1. The molten-pool problem. 



The portion of the above problem which is our focus is the conduction of heat 
in the solid, governed by the heat equation (Equation II. 2), to which the following 
boundary conditions apply: 



surface z = 0 


dT 

; k-^ = -q(r) 


(II.6) 


axis ;=0 


dr n 

: dr ~ 0 


(11.7) 


away /% z — y oo 


: T^T C , 


(11.8) 



where k is the thermal conductivity, and q(r) is the imposed surface heat flux (large 
at r = 0, falling off to zero at some small value of r, such that / 0 °° q{r)2irrdr = Q); the 
condition in Equation II. 7 results from the axisymmetric nature of the problem. Since 
the ultimate goal is to solve the Navier-Stokes equation in cylindrical coordinates, the 
heat equation will also be solved in this coordinate system. (see [Ref. 4]) 
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B. POINT-SOURCE PROBLEM 



Our goal is to Lake the first step toward solving the Xavier-Stokes equation 
( liquation 11.1) that arises from considering the molten-pool problem by solving the 
associated heat equation (Equation 11.2). We therefore assume cylindrical geome- 
try and azimuthal symmetry. We further simplify the conditions imposed on Equa- 
tion 11.2 by considering that the heat source applied to the horizontal surface of the 
block is a point source, with the total heat tlux Q — 1 . and that far away the solid 
approaches the uniform cold temperature of T c — 0. That is, the imposed heat flux 
'/(/*) = ^(r) is the Dirac delta function, with f£° r/(r) 2 rr 7 Y/r = 1. 

The existence of the solution to this type of problem is well established. How- 
ever. as is the case in general when Neumann boundary conditions apply, absent a 
specified initial temperature distribution, the solution is not unique (see [Ref. 5]). 
While we are interested in solving the point-source problem numerically in cylindrical 
coordinates, it has spherical symmetry, so it can be solved analytically in spherical 
coordinates more easily. Therefore, for the analytic solution, we rewrite Equation II. 2 
in spherically symmetric coordinates: 



8T 

dt 



K R?dR 




(11*9) 



where R is in spherical coordinates, R = \J r 2 + z l . Beginning with an initial tem- 
perature distribution of T = 0 everywhere and Q = 1, the solution to (II. 9) with 
boundary conditions 



is found to be 



surface z = 0 

axis r = 0 
far away r,z oo 

T(R, t) = — 



kj- = -<7M = S(r) 


(11.10) 


s-« 


(II.ll) 


T -* 0, 


(11.12) 


trfc Gv/d ’ 


(11.13) 
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where erfc(x) is the complementary error function. 

erfc(x) = l — erf(x), erf(x) = — [ e~ s ‘ 

x/ TT Jo 



els. 



As t — r oo, the steady solution becomes 



T(R) = 5 



] 



2ttkR 



(H.14) 



(for a derivation of this result, see [Ref. 5]). In (r, ~) coordinates, where R — v?' 2 + ~ 2 . 
these become 



T(r,z,t) = 



T(r,z) = 



2ttk\/ i’ 2 + ~ 2 



=erfc 



v/r 2 + z 
, '2\Zhct 



or 



(11.15) 

(11.16) 



2itk\/ r 2 + 

For computational purposes, the idealized problem of a semi-infinite solid 
is truncated to a finite domain in cylindrical coordinates with azimuthal symme- 
try. At the far boundaries on this finite domain, homogeneous Dirichlet conditions 
are imposed to approximate the conditions in the unbounded problem. To further 
simplify computational work, we assume that the (r, z) domain is the unit grid 
il = [0, l]x[0, 1] € 7£ 2 ;we will eventually assume k = 1. Thus, the equation we 
wish to solve is 



8T 

dt 



= kV 2 T , 



(11.17) 



with boundary conditions 



surface z — 0 : 


dT 


(11.18) 


axis r = 0 


: 

or 


(11.19) 


far boundaries v,z = 1 


: T = 0. 


(11.20) 
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III. 



DISCRETIZATION 



In order to solve the conduction problem numerically (Equation 11.2. with 
boundary conditions Equations IMS. 11.19. and 11.20). the equation representing the 
continuum problem must be discretized. That is. the values of the unknown function 
F art* to be determined from a set of equations which in some sense approximate 
Equation 11.2. Ultimately, a discrete approximation to the heat equation should meet 
t he following criteria: 

1 ) Provide for a unique solution to the problem. 

2) The solution should be '’( lose*’ to the exact solution for ‘‘sufficiently small" 
grid spacings. 

3) The solution should be "effectively computable." 

The significance of property 1) is obvious; property 2) relates to the questions of 
convergence and consistency for the discretization scheme. Property 3) relates to: 
a) the amount of computational work required to solve the problem, and b) the 
behavior of roundoff errors in the computed solution. The growth of the roundoff 
error is related to the notion of the stability of the discretization scheme, (see [Ref. 

«]) 

Solution of the heat equation (Equation II. 2) requires treatment of both space 
and time derivatives. The Finite Volume Element (FVE) method is used to discretize 
the spatial portion of the problem; finite differences are used to discretize the tem- 
poral portion. The approximation properties of the FVE method are fundamentally 
different from those associated with the finite difference method. Hence, this ap- 
proach treats time and space differently. The goal in applying the FVE method to 
produce spatial discretization is to make use of the advantages of the Finite Volume 
(FV) method, its ability to be faithful to the physics in general and conservation 
in particular, coupled with the flexibility of the finite element representation of the 
unknown functions (see [Ref. 1]). As outlined in [Ref. 1], the basic approach is 



to partition the spatial grid in two ways: as the union of a set of finite elements, 
the vertices of which comprise the grid points on which the unknowns are defined 
(see Figure 2), and as the union of a set of control volumes (see Figure 2), one for 
each grid point (the boundary points require separate treatment, as will be discussed 
later ). The elements are used to form basis functions, a linear combination of which is 
used to approximate the unknown function. Upon substitution of the approximation 
into i he equation to be solved, integrals over each control volume are taken. The 
integrals enforce conservation on each control volume, and therefore the partition 
enforces conservation on the entire domain. The system of equations generated by 
integration yields the discretization of the problem. 




0 1 2 N 

r ^ 



Figure 2. Partitioning of problem domain cross-section fi by elements. 

The process of discretizing Equation II. 2 in both space and time is complicated, 
the application of the FVE method to the space derivatives being the more complex 
portion of this process. Therefore, before proceeding to a discussion of the FVE 
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r ^ 

Figure 3. Partitioning of problem domain cross-section f2 into sub-volumes. 



method, the finite difference method is applied to the time derivative. In order to 
facilitate this presentation, we present a one-dimensional example of this type of 
discretization (see [Ref. 7]). Consider the one-dimensional version of Equation II. 2. 

c)T 0 2 T 



dt dx 2 ' 

where u = 1. One finite-difference approximation to Equation III. 1 is 

Tk, n+ l — Tk.n _ Tk+Un — 2Tfc,„ + Tfc_i,„ 

9 V 

where T is the exact solution of the approximating difference equations, 

H. = kh and t n = gn , (A;, n = 0, 1, 2, • • •). 

Let r = A, and this can be written as 



7fc,n+ 1 = l'Tk-l,n — (1 — '-l')Tk,n + vTk+\, n - 



(III.l) 
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Thus, we ('an compute the unknown 7V n+1 at the (n + l)st time step using the 
known values from the nth time step, giving an explicit formula for determing the 
unknowns. [Ref. 7] 

While this explicit relation provides a simple means to compute unknown 
values, is has a serious drawback. The process is only valid for 0 < & < i or 
// < . and therefore the time step Ai = g is necessarily small. Moreover, in order 

to achieve reasonable accuracy. Ax = h must also be kept small. There is, however. 

method which reduces the total amount of calculation and is valid (i.e., convergent 
and stable) for all finite values of /*, which was proposed by Crank and Nicholson 
(see [Ref. 7j). The technique is to consider Equation III. 1 as being satisfied at 
the midpoint {kh.(n + - 2 )g} and replace -r-jj by the average of its finite-difference 
approximations at the nth and (n + l)st time steps. That is, the equation 

dT\ _ (d 2 T\ 
dt /*,„+! “ 

is approximated by 

T^n+l — Tkji _ 1_ ( Tk-Un+l — 2Tfc, n +i + Tk+Un+l Tk-Un ~ 2Tfc iK + Tfc+I.n 
g 2 1 h 2 h 2 

which gives 

— r Tk-\,n+i + (2 + 2r)Tk <n +i — rT k+l , n+ i = 1 'Tk- t, n + (2 — 2 r)T k ,n + r T k +i,ni 

where r = -fa. Therefore, we now have three unknowns at the (n + l)st time step 
written in terms of three known values at the nth time step. Thus, the Crank- 
Nicholson method generates a set of equations which must be solved simultaneously; 
it is an implicit method. 

Although the Crank-Nicholson method for Equation III.l is stable for all pos- 
itive values of r in the sense that all errors eventually tend to zero as n tends to 
infinity, large values of r, e.g., 40, can introduce oscillations in the solution (see [Ref. 
7]). Therefore, a more general finite-difference approximation to Equation III.l is 
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tumid by using a weighted average of the finite-difference approximations of r— 7 at 
the //lii and (// + 1 1st time steps, which is given by 



I k.n + 1 F t . 71 



= n 



T, 



k-\ .iH-l 



~ 2Tk %n + 1 + 7fc+l. 



714-1 



/V 



+ (•- 



a 



Tk-l, n ~ -Tit, n T 7\*+ Ul 



where l) < a < 1 is possible. Note that o = 0 gives the explicit scheme, o — * , Crank- 
Nieholson. and a = 1 a fully implicit backward time-difference method. The scheme 
is unconditionally valid (stable and convergent) for 7 < a < 1, but for 0 < a < - we 
must have r < 1/(2 — 4a). Thus, for example, for a = 0, r — or (j < ~ , is a 

requirement to guarantee validitv.[Ref. 7] 

We now apply this finite-difference method of discretizing the time derivative 
to Equation 11.2. Following the work done above, a general finite-difference approxi- 
mation for the time derivative of 



ST 

dt 



kV 2 T 



(HI*2) 



is given by 



'T'n-H __ rpri 

= q«V 2 T” +i + (1 - a)KV 2 T n , 
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(HI-3) 



where g is the time step, and 0 < tv < 1 allows for a weighted average of the current 
and future time steps. The current time step is designated by n, n + 1 is the next 
time step; the value of a determines the type of time stepping. However, we have 
determined that in order to guarantee the validity (stability and convergence) of the 
explicit version of this scheme (a = 0), we must have g < We consider this 
value of g to be too small to be of computational use (which is verified by numerical 
experiments), and therefore we will use the values a = | and a = 1. Thus, we have 
the semi-discrete relationship between the current and subsequent time steps: 

(- - a«V 2 )T n+1 = (- + (1 - a)nV 2 )T n . (III.4) 

9 9 



This relationship is used to establish a discrete set of equations, which gives rise to 
the discretization of the problem. 
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A. BASIS FUNCTIONS: FINITE ELEMENTS 

\YV begin the discussion of discretizing Equation 1 1.2 in the spatial variables by 
recognizing that it is necessary to make an assumption about what types of functions 
may be used to approximate the unknown function T: we consider that T can be 
approximated by continuous, piecewise linear functions. A basis is then chosen for 
the space in which these functions are found. That is, a set of functions is selected 
whose linear combinations determine the space of functions of interest. Therefore, we 
choose basis functions. r, c), which are assumed to be continuous and piecewise 
linear, with 

r(r, -- 1 = ££;. ('••■)■ (III.5) 

k = 0 j= 0 

The next step is to determine how to construct these basis functions; the element 
partition of the domain 12 is used as a framework for this purpose, (see Figure 2). 

Since we are assuming cylindrical geometry with azimuthal symmetry, the 
directions of interest are r and a uniform grid of step size h = -jj ( N = number 
of grid points) is applied in both directions. In this way, we can designate function 
values on the grid 12 by 

T k ,] = T(rj,z k ), 

where = kh, Vj = jh , and k,j are the row and column indices. Each square 
is subdivided into two triangular elements by a diagonal (oriented in the direction 
of increasing r + z) 1 ; the basis functions to be used in approximating the unknown 
function are constructed using these triangular elements. For each interior grid point, 
or node 2 , there are six associated triangular elements: NNE, ENE, SE, SSW,WSW, 
NW (see Figure 4). 

'The orientation of the elements is based on the necessity to apply this technique to solving the 
molten-pool problem. In particular, we will be required to track the movement of the phase-change 
boundary between solid and liquid; element orientation has been chosen to facilitate keeping track 
of this moving boundary. 

2 For a two-dimensional grid with N intervals in each dimension, there are (N — l) 2 interior grid 
points, (N + l) 2 total nodes. 
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fk+lj+l) 



'M.) 




Figure !. Hat function, top view, k is the row index (z direction), j is the column 
index (r direction), y> direction is into the plane of the paper (0). 



If we suppose that this collection of elements is joined along shared borders, 
and that the elements can be “stretched’, then the Unite element for a particular 
node, say (A;, j ) , is formed by anchoring the collection along its external boundary (in 
the (/', z) plane) and extruding the grid point in a direction perpendicular to the (r, z) 
plain* and parallel to the tangent to the <p direction at the grid point ( k.j ). forming 
the familiar "hat" function (see Figure 5). Note that 



= < 



1 at the node (/•_,, Z*..) 
0 at all other nodes 



and is piecewise linear over each triangular element, giving the hat shape. 

By choosing these hat functions to be our basis functions, the coefficients of 
the dk,j (r, z) in Equation III. 5 become the values /3k, j = Tk,j. That is, we can now 
specify the nodal values of T on the grid as = T(r ; , z*,), with 



nv)-EEU^r4 

k=0 j=0 



(HI-6) 
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Since the Qk,j are piecewise linear, a continuous representation for function values of 
T is found by linear interpolation of nodal values between nodes. For example, 

Tk,j + Tkj+i 

2 ' 

Up to now, we have considered the sampled values of the unknown T at the 
(A r 4- l) 2 grid points as an (A r + l)x(N + 1) two-dimensional array; it will later 
be convenient to consider T as a one dimensional array. One simple method to 
accomplish this is to order the elements lexicographically. For example, 





© o 


To,, ••• 

r ltl • • • 


To,n+i ^ 
T\,n+i 






ran be written as 


V Tn+ l.o 


T n +i,i ••• 


Tn+ 1,N+1 j 






(r 0f o, 2o,i, . . . , Tojv+i, Ti t o, Ti,i, 


• • • ) T’l.N+l, • 


■ • , Tn+ 1 , 0 , rjv+1,1, • 


• • , ryv + i,yv+i) T , 
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which we will normally consider to he a column vector. Thus, we can relabel the 
vnlues T . ,i = /'/ . / = ( ,V + 1 )k -4- j + 1 and A :.j = 0 : .V + 1. in Equation III. 6 so that 

(/v+n 2 

T(r.z) = Ti<i>i{>\ z). fill. 7) 

;=i 

B. FVE STENCILS 

With the construction of the basis functions over the finite elements, we can 
now incorporate the control volumes to complete the finite volume element discretiza- 
tion. Below is an example the FVE method applied to a simpler problem: the purpose 
is to motivate and explain the work that is necessary to handle the specifics of apply- 
ing the FVE method to Equation 11.2. This explanation closely follows the example 
in [Ref. 1], 

1. Example: 2— D Potential Flow 

The example problem we consider is the potential flow problem in Euclidean 
two-dimensional coordinates. The domain for the problem is the unit grid $2 = 
[0. l]x[0, 1] € TZ i , where a: = 0 and y = 0 are the near boundaries and x = 1 and 
y = 1 are the far boundaries. The governing equation is 

V • (pWO = 77, (III.S) 

where p is the density, rp is the flow potential (i.e.. 0 X and xp y are the components 
of the fluid velocity in the x and y directions), and ij is the interior source flow rate. 
The boundary conditions are 

near boundaries : (pV?/>) • fi = Vh (Neumann boundary condition), 

far boundaries : rp = rpo (Dirichlet boundary condition), 

where h is the outward normal. Suppose that V is one of the control volumes in ft. 
similar to those depicted in Figure 3, where r and z are replaced by x and y. Since 
this problem is in two dimensions, the control volumes V are actually areas, and the 
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rurrospondine; surfaces S are the perimeters of these areas. Integrating Equation III. 8 
over l , we liave 

[ V-{pV*l>)dV= [ r] dV. 

Jv Jv 

By applying the (lauss Divergence Theorem, the volume integral on the left-hand 
side becomes a surface integral. 




(III.9) 



In terms of the physical units associated with the potential flow, the relationship is 



mass length mass 

— - — : area — — - volume; 

volume time time • volume 



each side represents a flow rate of mass per unit time, and (pVib) • h represents a flux 
across S. Therefore, the net flow from the interior source in V is balanced by the 
flow across the surface S, which can be interpreted as a mass conservation law for the 
control volume VC [Ref. 1] 

By partitioning the domain fl into a finite collection of control volumes (see 
Figure fi, where r and z are replaced by x and y), we can impose on each control 
volume the condition represented in Equation III. 9 (the boundaries require special 
treatment; a discussion of these considerations is included later). As noted earlier, 
the imposition of conservation on each control volume results in conservation over 
the entire domain. This integral condition on each control volume will be used to 
compute the discretization; the number of discrete equations is the number of control 
volumes, say m, used to partition f 7. To complete the discretization process, the 
unknown function ^ will be approximated by (3 in TZ m . By appropriately choosing 
basis functions, we can construct an approximation v expressed in terms of its nodal 
values. Consider a triangular element partition similar to that in Figure 2 (where r 
and r are replaced by x and y ) and let T be the space of continuous, piecewise linear 
functions on fl associated with this partition. Ignoring for the moment conditions at 
the boundaries, the FVE discretization of Equation III. 8 comes from finding a v € T 
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-ii( Ii that lunation 1 1 1.9 holds for uadi of the control volumes in the partition of fh 
1 lie problem in R m is then defined when e is expressed in terms of nodal basis for T: 



— y ^ 



HI. 10) 



w=\ 



where u w is the value of v at the w th node and o w is the “hat function ' associated with 
the wth node (as above, we lexicographically order the nodes, resulting in a single 
subscript). Substituting Equation III. 10 into Equation III. 9, we have the matrix 
equal ion 



Lu = /, 



where L is the nxn matrix with entries 



L, u ,z = / {pV6 z ) 



i dS 



(III. 



( 111 . 12 ) 



ami 



f w = [ vdv (HI. 13) 

(except for boundary conditions). [Ref. 1] 

The treatment of the boundaries depends the type of boundary condition. 
Neumann conditions can be imposed indirectly by substituting the flux value V’i into 
the appropriate term in Equation III. 9. Specifically, suppose we choose the control 
volume V that includes the origin (lower left-hand corner, Figure 3), where the surface 
of the control volume coincides with the boundary along the two axes. The integral 
condition for V is 



[ fadS + / (pVxb) • hdS = / r/dV; 

Jx y y = 0 Jx y y^0 Jv 



(III. 14) 



Equation III. 14 is substituted for the interior condition in Equation III. 9 into the 
discrete approximation v for the corner volume containing the origin. [Ref. 1] 

Dirichlet conditions are imposed directly; u w takes on the value of t/q at each 
Dirichlet node, i.e., each node on the far boundaries (x,y = 1). This approach may 
lead to fewer unknowns u w than equations, a problem easily resolved by discarding the 
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('({nations associated with the Dirichlet nodes. In this case, the collection of control 
volumes no longer partitions the domain and. therefore, conservation no longer is 
strictly enforced. However, it is at the Dirichlet nodes where the loss of conservation 
is not a concern in general. [Ref. 1] 

Assuming that there is a balance in the number of equations and unknowns, we 
now must implement a numerical rule for evaluating the integrals in Equations III. 12 
and III. Id. For Equation III. 12. we use the following rule on each segment A that is 
part of the interior surface S associated with a volume V: 



[ (pVd>) • hdS ~ p( ) [ V&'hdA 
J a J A 



i a J A 

where P 4 is point of intersection of A and the grid lines passing through the node 
associated with V. For Neumann boundary segments, we use the quadrature rule 

[ fadS ~ v>\(Pa)\A\, 

J A 

where |. 4 | is the “surface area”, i.e., length, of A. For Equation III. 13 , we use 

/ r,dV~ii(N v ) \V\, 

J V 

where Nv is the node associated with V and |V| = f v dV is the “volume”, i.e., area, 
of V\[Ref. 1] 

Numerically evaluating Equations III. 12 and III . 13 yields the FVE discretiza- 
tion of Equation III. 8 . The value associated with each node, say (p,q), is a sum 



w=p+\ ,z=q + l 



L 



&W,2 W,Z ? 



w—p— l ,z=q — I 

where the coefficients a Wy2 are determined by the integration. A convenient way to 
represent the sum associated with each grid point is to place the coefficients a Wi2 into 
a stencil, 



a p+ 1,<7 

C*p,q— I ^p,g ^p,<?+l 
I,<? 



^P»<? — ®p— 1,9 u p— 1 ,q d~ — 1 U p>q— l 

d" ^p,<7^p,<? d" <^p+l ,<?^p+l,7 d“ C*p,g + l W p,<? + 1 
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where. fur this example, the corner values in the stencil are all zero. The value for 
iudr \ !>.({) is determined by applying the stencil to u ; , 7 . 

To see what type of stencils are produced, consider the discretization on a 
uniform mesh with grid size ^ ~ in both coordinates. We use double subscripts 

/).</ — 0 : m — 1, where 0 corresponds to the Neumann boundary nodes and rn — 1 
corresponds to the Dirichlet boundary nodes. Written in stencil form, the interior 
nodes. 0 < p,q < rn — l. for the equations in Equation III. 11 are as follows ( Y 2 
indicates the sum of the outer stencil entries): 



p{{p + r 2 )h,q ll) 

p{pi> ■ (<i - y/o -e f>(p h - (<i + i)/o 

p{(p - 

For the corner Neumann node (p, q = 0), the stencil is given by 



u p q - k 2 ij{pli. (ih). 



ip(io) 

-E *p(o*) 



h 1 x h ( , , h. h \ 

uo.o = — r;(0,0) - - 0i(O, -) + t/>t(-,0) . 



For both stencils, the coefficient of the unknown to which the stencil is applied is 

-E- 



For grid points adjacent to a Dirichlet node, the stencil entry reaching to the 
Dirichlet node value is moved to the right-hand side as a coefficient of the boundary 
data. Otherwise, the equations for these points are the same as above. For example, 
at the point (0 < p < rn — 1, q = m — 1) we have 



p{{p + \)li, 1 - h) 

-E 

p((p - |)M - h ) 



Up,tn- 1 = h' 2 7](pll. 1 - ll) 



-p(ph, 1 - (ph, 1), 



where is the sum of the outer stencil entries without the boundary terms removed, 

= F(A 1 ~ + P(ph, {q ~ |)/») + p((p + ^ )h , 1 -h)+ p{{p - ^)h, 1 - h). 
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(see | Ref. 1]) 

2. Application to Axisymmetric Heat Equation 

This technique is now used to compute the FVE stencils for Equation 1 1.2, bv 
combining the finite elements (see Figure 2) with the control volumes (see Figure 3). 
A control volume for the Zth grid point, call it IT is a toroidal prism (see Figure 6). 
It. is generated by taking the two-dimensional sub-volume for that point (in the (r, r) 
plane), and sweeping through 360° in the p direction. 

♦ z 



y 

x r 




control 

volume 



Figure 6. Toroidal volume for the conduction problem. 



Integrating Equation III. 4 over all control volumes Vj, where V = Ui 
partitions the domain fi, we have 



(AM- 1) 2 



Vi 



T Q - a/cV 2 ^ T n+l dV = T Q + (1 - a)fcV 2 j T n dV , (III-15) 
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or, upon application of the (iauss Divergence Theorem the volume integral of the 
V 2 T term becomes a surface integral, and we have 



T n +'<iv - aK [ vr n+1 . hds = - [ r n dv + (i - <v)« [ vr n • hds. (iii.ig) 

(J JV Js (J Jv Js 

where h is the outward normal. Substituting the expression for approximating the 



unknown function T (Equation III. 7). we have 

- [ <j>? +l dV - an [ V<ft +l • fids 

(I Jv Js 



(;V+I > 2 r 

E 



/=I 



is 

(iV+l ) 2 



T" +I 



(III.17) 



= E 

/= 1 



- [ <j>?dV + (1 - «)k [ Vd>” ■ hdS 

a Jv Js 



T 'Tl 

I ’ 



1 

J) JV JS 

where we now have integrals over each of the (A + 1)"’ control volumes, resulting in 
a set of ( A’ + l) 2 equations. This set of equations can be written both as an operator 
equation. L[T n+l ] = /(T n ), and as a matrix equation 3 , M(T’ n+l ] = /(T n ), where the 
operator L and the matrix M are given by 

'1 



L 

M 






= / i- 

Jv \9 

= - [ tf +l dV -qk f V4>? +i ■ hdS , 
g Jv Js 

f{T n ) and f{T n ) are given by 

/(T n ) = J v Q + (i -«)KV 2 jrw, 

(N+1)2 n r r 

f{T n ) = E -/ tfdV + (1 - <x)k I yw-h 

I T Q JV J S 



77. 



Any function whose coefficients satisfy the resulting set of discrete equations neces- 
sarily satisfies the conservation law over any volume made up of the union of control 
volumes, except possibly at the boundaries. The Neumann conditions at r = 0 and 
z = 0 are incorporated indirectly into T, by substituting the (zero) normal derivative 



3 The operator may be represented as a continuous operator, L; as a discrete operator, L h \ or as 
a matrix, M. 



21 



value into the appropriate term in Equation III. 15. The Dirichlet conditions at the 
lar boundaries are imposed directly on T: control volumes that would usually be as- 
sociated with these boundary nodes are eliminated (see Figure 7). Thus, the reduced 
collection of control volumes no longer partitions the grid fl. which slightly impairs 
conservation in the discretization. However, since we require the temperature to fall 
oil to zero at these points, this loss of conservation is not a concern. [Ref. 1] 




0 1 2 N 

r 



Figure 7. The problem domain cross-section is depicted after partitioning into sub- 
volumes. Homogeneous Dirichlet boundary conditions obviate the necessity to define 
volumes for grid points at the far boundaries. 

The integrals over the basis functions have been computed using Maple soft- 
ware and a program designed by David Canright (see Appendix B and [Ref. 8]). In 
order to calculate the integrals, we must be able to represent the unknown function 
over the basis functions, which themselves are collections of triangular planes. There- 
fore, the method is to determine first the function for the plane through three given 
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points (the rases of vertical planes and three points collinear are not considered). 

I his procedure is then used to create the triangular elements which, when assem- 
bled. form the hat function (see Figure 5). For a given grid point, sav (/j, j), six of 
these triangular planes are joined, corresponding to the six triangles surrounding the 
point. XNE. FNE, SE. SSW.WSW, NW (see Figure 4). The unknown function is 
then interpolated over the six elements. 

Once the unknown function is represented by the combination of these six 
interpolations, we can use Maple to compute the volume integrals, and the surface 
integrals of the gradients, in Equation III. 17. The results of these integrals provide 
the coefficients which comprise the FVE stencils. 

rhus, we have the following stencils for the volume and surface integrals re- 
spectively for interior grid points: 



L 






(Z~£T r .,^,,)dv = ^ 



32 j — 5 16j + 5 

32 j — ll 224 j 32j + 1 1 

16 j — 5 32 j + 5 



T, 






(III. 18) 



and 



JS *.J p q Z 



2 j 

2] - 1 - 8 j 2j + 1 



Tk,j, 



(III. 19) 



2i 



where the 2n resulting from integration in the < p direction has been factored out. 
(Recall that in cylindrical coordinates, 



dV = j U f J+i f Zk * h rdzdrdy. 

Jo J r 1 J z, i 
k ~i 

Note that control volumes increase with radial distance from the origin, which is 
reflected in a radial bias in the stencils. More specifically, rj = jh is the radial 
distance from the origin, with j the column index for the unknown matrix, h the 
meshsize. Thus, stencil values increase with distance from the origin. 

These stencils are applied to the unknown at a point on the grid, designated 
by its row/column position (Ar,j). Stencil entries indicate the values to be applied; 
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blank entries indicate a value of zero. Entry position indicates to which grid point 
i lie value is applied: left/right and up/down in the stencil correspond to neighboring 
points in those directions on the grid. That is, if stencil values are replaced by the 
positions to which they apply, we have 

(k + l.p) (k + \ ,j + 1) 

(k.j — I) (k'j) (k- J + 1) 

(k - \ .j - \ ) {k - \ .j) 



rims, fur example, the value that appears in the (k + 1, j) stencil position is applied 
to the grid point that occupies that same position. 

Using technicpies similar to those outlined in the two-dimensional example, 
boundary point stencils have been computed for the volume and surface integrals 
respectively. At the origin, r = 0, z = 0: 



_f_ 

384 





1 5 




l 




15 3 


J 

7o,o, and — 

O 


-3 2 












(III. 20) 



at r — 0: 



384 



1 


5 




and — 

8 


1 


25 


11 


Tk,o, 


-6 4 


6 








1 



Tk, o; 



(III. 21) 



and at z — 0: 



If 

384 



32 j - 5 1 6j + 5 

24 j — 8 1 12 j -(-5 8j T 3 



Toj, and — 



4j 

2j - 1 -8 j 2 j + 



To j. 



{ 11122 ) 

For points adjacent to the far boundaries, the stencils in Equations III. 18, III. 19, III. 21, 
and 1 1 1.22 are applied. However, since the homogeneous Dirichlet conditions dictate 
that these boundary values are zero, the resulting contribution after stencil appli- 
cation remains zero. Thus, in effect, far-boundary values do not contribute to the 
stencils for points adjacent to these boundaries. 
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We ran now combine all of the above stencils to generate the operator L h , 



L k [T n+ '} = f k (T n ), 



(1 11.23) 



where li is the step size on the grid. The matrix representation of L h . M. is a block 
tridiagonal (.V + 1 ) ’ x ( A'’ + l) 2 matrix of the form 



M = 



A A 0 0 

r .4 A 0 

0 r .4 A 

0 0 •• 



0 



0 0 
0 0 
0 0 
0 



0 0 r .4 A 

o o r .4 



where .4, A, and T are (N + l)x(lV + 1) generic banded matrices (not all identical); .4 
is tridiagonal. A is upper bidiagonal, and T is lower bidiagonal. When M multiplies 
the matrix of (N + l) 2 values of T, arranged lexicographically as a column vector, 
the matrices A , A, and T produce the effect of the stencils “reaching’’ function val- 
ues respectively on the current row, the row above, and the row below. Numerical 
experiments using Matlab to construct the matrix M for N = 8,12,16,20, and 24, 
with g = li and a = and 1, indicate that it has full rank. Thus, we expect a unique 
solution to the linear algebra problem which is a discretization of the heat equation. 
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IV. 



RELAXATION SCHEMES 



T lie FYE method, with weighted-average time-stepping, has 


been used to 


discretize the continuum problem. 




-N 

[> 

II 

f- | 


(IV.! ) 


<)t 


by 




L h [T n+ 1] = f h (T n ) 


( IV. 2) 


or. in matrix form. 




M[r ,+l ] = /(f n ), 


(IV.;n 



on a grid of meshsize It — The input for these equations is the solution at the 
current time step, T n (where T n and T n are used interchangeably), the result of 
solving the equations is the solution at the next time step, T 714 " 1 . As was indicated 
in Chapter III, there is good reason to believe that the matrix M is of full rank and, 
thus, expect a unique solution to exist for the linear algebra problem that arises from 
the discretization of the continuum problem. Solution by direct methods requires 
factorization of M and, since M is both large and sparse, solution by direct methods 
may be impractical. We therefore consider iterative methods to solve the matrix 
equation at each time step. These methods generate, for each time step, a sequence 
of approximate solutions, {T^j} 1 ; the choice of relaxation scheme determines whether 
or not the sequence {T( 5 )} converges. Upon implementation of a relaxation scheme 
that produces a converging sequence of approximate solutions at each time step, we 
begin with an initial temperature distribution, T°, and the solutions are stepped in 
time. The desired result is a sequence of solutions, corresponding to a sequence of 
t ime steps, which tends to steady state. This process allows for the evaluation of 

1 Here the subscript (.9) indexes the sequence of approximate solutions; elsewhere, subscripts 
without parentheses indicate grid position or vector components. 
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various values for the time step, as well as of various weightings used in averaging the 
current and subsequent time steps. 

A. ITERATION MATRICES 

It is common in constructing iterative methods to propose a splitting for the 
matrix M = E — F. where linear systems of the form Ex = b are “easy" to solve (see 
Kef. !)]). The sequence of approximations. {T (s )}, is generated hv ET( S+1 ) = FT (s ) + / 
and. therefore, it is natural to define an iteration matrix. P = E“'F. so that 
f ls+ i) = P T {s] + E 1 /. Additionally, if T* is the exact solution so that MT* = /, 
then E T* = FT* + / and T* = E _1 FT* + E ~ l f. Hence the solution. T*. is a fixed 
point of the iteration T( s+1 j = E -1 FT( S ) + E ~ l f. 

The matrix P = E -1 F is also called the error propagation matrix, since if 
T (s+ |) = PT( s j + E -1 / or T( J+1 ) = PT( S ) + B (B a constant vector), and e( s+1 ) is the 
error at the (s + l)st step, then 

e(s+i) = T( s+1 ) — T* 

= [P f (s) + B\ - [Pf * + B) (since f* = PT* + B) 

= Pf (s) -PT* 

= P (f (s) -f*), 

and therefore 

e (s+ i) = Pe (s) . (IV.4) 

Using induction, it is easy to show that e( s ) = P 3 e(o). By Equation IV.4, ep) = Pe(o) 
and 



e (2 ) = Pe ( i) 

= P[Pe ( o)] 

= P 2 e ( o)- 
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My induction, assume that t^) = P fc C(op in order to show that c ^ + n = P^'bfo). 
Now. 



‘V+ 1) = p c (A;) 

= P[PV (0) j (by the induction hypothesis) 

= P ( "' +I, "(o,. 

Thus, since e^+i) = P ,A, ’ + 'Ie(o), we conclude C( s ) = P S C( 0 ), for .s any integer. Addi- 



i ionallv. 



Ie“ II = | ! P S 



(«))H E 



< I'PI 



(o)l 



(IV. 5) 



which will he useful later. 

One of the simplest relaxation methods is the Jacobi method, also referred 
to as simultaneous displacement (for a more detailed account of all of these methods, 
see (Ref. 3] or [Ref. 9]). It is produced by solving the / th equation of Equation IV. 3 
for the /th unknown, T), / = 1 : (N + l) 2 . Before proceeding to a discussion of the 
iteration matrix, we present the component form for this iteration scheme: 



rp(new) -p(new ) 

1 kj ~ 1 l 



h - )i 

^224, + *±8j 



( IV .6) 



where Ylv anc ^ are * respectively, the sum of volume stencil entries and surface 
stencil entries applied to the current values of the unknown, except at the grid 
point on which the stencil is centered: 



/ 



old 



E, = 



m-yTg* +(i6j + 5 )r;;'“’ +1 



+(32j- 

+(16J - 5) TltZ-1 +(32 j + 5 )TltZ 



(old) \ 
■ 1.1 + 

+(32; + 11)7^ 

) 



(IV. 7) 



uid 



* old 

E, = 



2 jlitu 

+(-1+2 +(>+2 
+ 2jT<"» 



(IV.8) 
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OiK‘ iteration sweep consists of computing Equation IV. 6 for each component of the 
unknown vector T. Provided that the process converges, the sweeps continue until a 
desired level of convergence is reached. 

Another, more succinct, way to present this method is in matrix form. If 
we write the operator matrix as the sum of a diagonal matrix (D), a strictly upper 
triangular matrix (U), and a stictly lower triangular matrix (L), 

M = D + U + L, 

the matrix equation to be solved becomes 

(D + U + L)[fj=/. 

Now define Ej = D and F j = — (U + L) , so that this may be written as 

D[f] = -(U + L)(f] + /, or 
E j[f] = F4f] + /, 

or as 



T = -D-‘(U + L)[f] + D- 1 /, or 

f = Ej'Fjlfl + Ej'/, 

which corresponds to solving the l th equation for T}, / = 1 : (N + l) 2 . If we define 
the .Jacobi iteration matrix as 



P J — Ej’Fj, or 

= — D -1 (U + L), 

the matrix form of the Jacobi method becomes 

f(neur) _ p^f{old) ] + gjt J (W.9) 
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\ modified version of this method, called the weighted Jacobi method, is 
determined by introducing a weight. U £ ^ 1 (*j = 1 is the original .Jacobi method) 



where I is the indentity matrix. = D. and = ( I — lj)D — ~’(U + L), and 



is the weighted Jacobi iteration matrix. In this method, the new approximation is a 
weighted average of an intermediate approximation and the old approximation: the 
intermediate approximation is the Jacobi iterate of the old approximation. 

The weighted Jacobi method computes all of the components of the new ap- 
proximation before it begins to use them in the next approximation. This requires 
storing both the current and new approximations; a simple modification allows for 
updating the current approximation “in place", using only the storage required for 
one approximation. The Gauss-Seidel method incorporates the new changes as soon 
as they are computed by overwriting the current approximation with the new (see 
[Kef. .*}]). More important, however, is that (on model problems) the Gauss-Seidel 
method generally converges about twice as fast at the Jacobi method (see [Ref. 9]). 
In component form, this iteration scheme is 



(see I Kef. .$]). The matrix form is 



'p(new) p 




(IV. 10) 



P = E"'F 

1 u; — -A jJ 

= D“‘[(l — lj)D — u,’(U + L)] 

= ( 1 — cJ)I + tv'P J 




(IV.ll) 



where the arrow indicates overwriting, and 





+(32j + ll)7t“> , (IV. 12) 
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nut 






V" : 
z— s 



+(-l +->j)T^ 



neiu) 
1 






-J 1 k + 1 j 



, •> -rp(new) 

+ ~J 1 k-\ ,j 



+(i j)Ti:‘i\ 



(IV. 1:!) 



/ 



III matrix form, using M = D 4 - U + L. Er; = (D + L), and Fr; = — U, we have 

(D + L)[f] = -U [f] + /, or 
E G [f] = F G [f\ + f, 



f(nrw) ^ _( D + L)- x U[f {old) ] + (D + L)- 1 /, or 

f(new) Ec l Fc[f {old) } + Ec l f. 



This corresponds to solving the l th equation for T), using the new approximations for 
components 1,2 — 1. Now, define the Gauss-Seidel iteration matrix, 



Pa = Ea'Fa 

= — (D + L)~'U, 



so that the iteration scheme in matrix form becomes 

f(new) p G [f(»W)] + E" 1 /. (IV. 14 ) 

With the above lexicographic ordering of the (N + l) 2 components of T and the 
components updated in ascending order, the effect of a Gauss-Seidel sweep is to start 
at r = 0 and update in the radial direction for each vertical step, starting at z = 0 
(see Figure 8). 

As with the weighted Jacobi method, we can make a modification to the Gauss- 
Seidel method. Define a parameter 7 £ 7 Z (7 = 1 is the original Gauss-Seidel method) 
and the modified method is successive over-relaxation, SOR, (see (Ref. 9 ]) which, 
in matrix form, is given by 

f(ner.) _ p ^f(otd) + ^-1 £ (IV. 15) 
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Figure 8. Gauss-Seidel sweep. 

where, with E 7 = D + 7L and F 7 = (1 — 7 )D — 7U, we have 

P, = (D + 1 L)- , [(l - 7 )D-tU] 

= 

Similar to the weighted Jacobi method, SOR is a weighted average of an intermediate 
approximation and the old approximation. 

The question of interest now is whether or not the sequence of approximations. 
{T( a .)}, generated by ET( S+1 ) = FT( S ) + /, converges. Therefore, we make use of the 
following theorem: 

Theorem IV. 1 Suppose that / £ 7 v n and M = E-F 6 7 Z nxn is nonsingular. If 
E is nonsingular and the spectral radius of E _1 F satisfies ^(E^F) < 1, then the 
iterates T( s p defined by ET( S+I j = FT( S ) + f , converge to T — M ~ l f for any starting 
vector T( o)* 

The proof is found in [Ref. 9], and makes use of the following lemma: 
Lemma IV. 1 If p( E _1 F) < I, then (E^F)* -4 0 as k — 00. 
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N = 8 | 


X = 12 


.V = 16 


X = 20 


X = 24 


•Jacobi 


.9259 1 


| .9544 | .9675 


.9749 


.9796 


cJ = 


.9721 


.9826 


.9875 


.9903 


.9921 


uj = - 

3 


.9442 


.9652 


.9750 


.9806 


.9841 


Gauss-Seidel 


.8395 


i .8982 


.9263 


.9425 


.9530 


7 = ^ 


.9670 


.9793 


.9851 


.9884 


.9905 


7 = 5 


.9189 


.9488 


.9630 


.9711 


.9764 



Table I. Spectral Radii for Crauk-Nicholson Time Stepping (a = ^). 
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X = 8 1 


t .Y = 12 ! .V = 16 


A r = 20 


A r = 24 


Jacobi 


.9499 i 


.9709 


.9801 


.9850 


.9881 


w = i 


.9817 


.9892 


.9925 


.9943 


.9955 


^ - 3 


.9634 


.9784 


.9850 


.9887 


.9909 


Gauss-Seidel 


.8931 


.9362 


.9556 


.9663 


.9730 


7 = ^ 


.9782 


.9871 


.9910 


.9932 


.9946 


7 = j 


.9462 


.9680 


.9777 


.9831 


.9865 



Table II. Spectral Radii for Implicit Time Stepping (a = 1 ). 

The matrices M, Ej, and Eg, and iteration matrices, P, for various values 
of iV, a = | and 1 , and g = /i, have been constructed using Matlab. We have 
experimentally verified that M, Ej, and Eg are nonsingular, and the spectral radii of 
the error propagation matrices, P. have been computed. Due to memory limitations 
with Matlab , the calculations are made for relatively small values of N. The results 
are indicated in Tables I and II, where the value of a indicates the type of time 
stepping, and the values of u and 7 are the weightings in the weighted Jacobi and 
SOR methods respectively. 

The results of these numerical experiments indicate that both p(Pj) < 1 and 
/?(Pg) < 1, with p(Pg) < p(Pj). The modifications to the Jacobi and Gauss-Seidel 
methods do not provide any improvement; the spectral radii for both methods are 
higher for the modified iteration matrices than for the original iteration matrices. 
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I 1ms. the ( lauss-Soidoi method, used with either the Crank-Nicholson or implicit 
i ime-si epping scheme. appears t o he the best choice of these relaxation schemes. How- 
ever. il the spectral radius of P is near unity, convergence may by unacceptably slow 
mice the error tends to 0 like/)(P)*’. as indicated bv the lemma and J|f 5 || _ P || s II (o) j| 
(liquation IV. 5). As is evident from the above tal)les. even for a moderately spaced 
grid (e.g., N = 16. (j = /?.. and o = \ or 1). p(P) > .9, and the spectral radius 
increases with A . 

B. ALTERNATIVE SCHEMES: LINE RELAXATION 

Both the .Jacobi and Gauss-Seidel methods give rise to iteration matrices, are 
implemented in a straightforward manner, and are attractive because of their sim- 
plicity. While the ease of implementation is a significant advantage, the convergence 
properties of either or both may not be acceptable and, therefore, alternative schemes 
may be desirable. One such, which seems reasonable based on the geometry of the 
problem, is line relaxation. There are two obvious options in this regard: radial or 
vertical line relaxation (see Figures 9 and 10), where either an entire row or column 
is updated simultaneously. Both options require the solution of an (A r + l)x(A r + 1) 
tridiagonal matrix for each row/c.olumn update in the unknown matrix. That is. while 
the previously outlined relaxation schemes (point relaxation) proceed by successively 
solving a collection of algebraic equations for one unknown, line relaxation proceeds 
by solving a succession of matrix equations. For example, suppose that the matrices 
.1, A, and F are tridiagonal, upper bidiagonal, and lower bidiagonal, respectively, and 
that T x and /,* are the portions of the respective unknown and right hand side vectors 
in Equation IV. 3 that are rows in their corresponding matrices: 
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Radial line relaxation consists of solving t lie following succession of matrix equations: 

T> A~ l {f\ — AT 2 ), 

- A-Hh-rr, - at 3 ), 

7s *- .A _, (/ 3 - VT 2 - AT 4 ), and 

<- .4-*(/4 - rr 3 ). 

In order to illustrate how line relaxation works relative to the FVE stencils, 



let 



row 
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(IV. 19) 
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-y> ( new fcoluj/m 
1 '.j 
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-i( nen 1 ) t 
A.' 4- l . j J 



m) t.lial. for the radial and vortical relaxations respectively, we must solve 

(IV. 20) 
( IV. 2 1 ) 

simultaneouslv for each row or column in the unknown matrix. 
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next row update 
latest row update 



r 



“ *“ — old approximation 
* new approximatton 

Figure 9. Updating an entire row at one time. 

The computational cost of updating an entire row/column at a time is higher 
than a row/column update using either the Jacobi or Gauss-Seidel methods. However, 
this type of relaxation may be sufficiently efficacious to warrant the extra expense, in 
that convergence may be achieved significantly faster. For comparison, we now con- 
sider local mode analysis of the Jacobi and Gauss-Seidel methods and line relaxation. 
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latest column update 



next column update 




r 

” ~ ” old approximation 
new approxunauon 

Figure 10. Updating an entire column at a time. 

C. LOCAL MODE ANALYSIS 

While an examination of the spectral radii of iteration (error progagation) 
matrices can be instructive, it is often useful to conduct a more detailed analysis. 
\\ hat follows is a local mode analysis of various iteration schemes. Since the error 
propagation matrices indicate how the error evolves during the iteration process, we 
use a DFT (see Appendix A, [Ref. 10]) to expand components of the error equations 
associated with the various relaxation schemes. The coefficients in these expansions 
are the factors by which the corresponding modes of the error are magnified/reduced 
for each relaxation sweep. Thus, by examining these coefficients, we can determine 
how quickly the error tends to zero, which indicates how quickly the solution con- 
verges. The following analysis does not apply to points on the boundary, it is only 
valid for interior points. Thus, while it gives more information about the functioning 
of the scheme in the interior of the domain, boundary peculiarities are not addressed. 
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We begin bv recalling that the current error. C( 5 ), is the difference between 
the exact solution. F*. and the current approximation. T( s ), 

f ( S ) T ^(s)* 



where we desire that the sequence. {c( s ) }, tend to zero. Substituting this expression 
into Equations IV.fi and IV. IK we have the following error equations for the Jacobi 
and Clauss-Seidel methods where, in order to avoid confusion with tin" exponential 
luncrion in the DFT expansion, we represent the components of e by r^ tJ : 



T i • (ueu/) 

Jacobi : = 



Gauss — Seidel : = 









IV. 22) 



(IV. 23) 



where Hv^ and VX are now applied to the (see Equations IV. 7. IV. 8. IV. 12 

and IV. 13). Equations IV. 22 and IV. 23 indicate how the sequence {e( s )} evolves 
during the iteration process. If we expand these relationships in a discrete Fourier 
transform, (DFT) (see Appendix A or [Ref. 10]), the magnitude of the transform co- 
efficients will indicate the “amount” of each mode of the error that is present in each 
component of ef( s ). We then compare coefficients to determine the ratio between the 
new and current approximations for each component of the error. In other words, the 
DFT allows us to analyze the growth of the error by examining component behavior. 

Expanding Equation IV. 22 in a DFT, where 
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C (/, rn) = e * e N 



^ 0, and F{j) - 



Sj224j + **8j' 



and the superscripts (o) and (n) indicate the components of the old and new approx- 
imations, we have 
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Making use of the orthogonality property of the complex exponential (see Appendix A) 
we <'an equate individual terms of the sum, and then divide by C(/,m) to give 
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N = i6 j 


N = 32 


1 


J = 1 


1 J = A i J = A' - 1 |( j = 1 


j 


J = v — 1 


l '- = 1 


1 .009 


1.009 ! 1.009 


I 1.026 


1 .026 


1 .026 
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0=1 
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1.013 


1.013 



Table III. Maximum Values of 
Relaxation. 



' for l, in = 0 : .V for the Jacobi Method of 
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In order to determine the greatest factor by which a mode of the error is multiplied. 



IT 



(")| 



we seek the maximum of ';™ ~ over the values /, m = — rr + 1 : ■+ . In other words, we 

IT, ml 

seek to determine a bound on how well we can expect this type of relaxation scheme 
to perform. This ratio is a function of N,l,m and j and, while difficult to determine 
analytically, may be calculated numerically. Using Matlab, the maximum of this ratio 
for /. m = 0 : /V has been calculated for A" = 16,32; a = y 1; and j = 1 , y , .V — 1 ; the 
results are indicated in Table III (see Appendix A for a discussion of the equivalence 



of centered indices, /,m = — — + 1 : y, and noil-centered indices, l,m — 0 : .V). 
Additionally, by considering the matrix of grid points /,m = 0 : N, we can determine 
a correspondence between sample points on the grid and type of associated frequency 
(see Figure 11 and Appendix A). The matrices of values of for j = N — 1 are 
depicted in Figures 12 and 13. 



The information in Table III indicates that the Jacobi method results in a 
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